Random variables

What you need to know

• Summing independent random variables results in another random variable called sumple sum. The mean of the sample sum is different from the population mean or expectation which is an exact quantity we want to approximate by sampling.

• The Law of Large Numbers is a principle that states that as the number \(N\), the sample mean approaches the population mean with a standard deviation falling off as \(N^{-1/2}\)

• The Central Limit Theorem (CLT) tells us that summing independent and identically distributed random variables with well-defined means and variances results in Gaussian distribution regardless of the nature of a random variable.

• A model of random walk describes the erratic, unpredictable motion of atoms and molecules, providing a fundamental model for diffusion processes and molecular motion in fluids. The number of steps to the right (or left) of a 1D random walker results in a binomial probability distribution. Following CLT binomial distribution in the large N limit can be shown to be well approximated by gaussian with the same mean and variance.

Introducing random variables

• A random variable X is a variable whose value depends on the realization of experiment or simulations.

  • \(X(\omega)\) is a function from possible outcomes of a sample space \(\omega \in \Omega\).

  • For a coin toss \(\Omega={H,T}\) \(X(H)=+1\) and \(X(T)=-1\). Every time the experiment is done, X returns either +1 or -1. We could also make functions of random variables, e.g., every time X=+1, we ear 25 cents, etc.

• Random variables are classified into two main types: discrete and continuous.

  • Discrete Random Variable: It assumes a number of distinct values. Discrete random variables are used to model scenarios where outcomes can be counted, such as the number of particles emitted by a radioactive source in a given time interval or the number of photons hitting a detector in a certain period.

  • Continuous Random Variable: It can take any value within a continuous range. These variables describe quantities that can vary smoothly, such as the position of a particle in space, the velocity of a molecule in a gas, or the energy levels of an atom.

Random numbers in python

• The numpy.random has the fastest random number generators based on low-level code written in C.

• The Scipy.stats has an extensive library of statistical distributions and tools for statistical analysis.

• First, we take a look at the most widely used random numbers of numpy, also called standard random numbers. These are rand (uniform random number on interval 0,1) and randn (stnadard average random number with 0 mean and 1 variance).

• When running code that uses random numbers results will always differ for every run. If you want code to reproduce the same result, you can fix the seed to get reproducible results: np.random.seed(8376743)

• To convert random variables to probability distributions we need to generate large enough sample then perform histogramming via np.hist or directly histogram and visualize by one shot via plt.hist()

import numpy as np
import matplotlib.pyplot as plt

X = np.random.rand(50)

print(X)
plt.plot(X, '-o')

[0.01649061 0.94809928 0.53167608 0.70297179 0.24938418 0.644101
 0.63191325 0.26702288 0.54860171 0.9113438  0.42463562 0.64045427
 0.5836673  0.87348577 0.87839506 0.30074767 0.99802114 0.8662603
 0.59231173 0.92788546 0.48898307 0.46747427 0.90064271 0.22878622
 0.50978301 0.18799548 0.80011062 0.19234139 0.37394562 0.75512087
 0.90381803 0.57959818 0.62468977 0.69036726 0.73657662 0.31476815
 0.51585703 0.55235361 0.03559257 0.95396299 0.75260844 0.9636808
 0.92556659 0.38376044 0.94806819 0.51243658 0.5512443  0.00192689
 0.69598612 0.21691083]

[<matplotlib.lines.Line2D at 0x7feac5bb1040>]

Probability Distribution of a Random Variable

• For any random variable \( X \), we are interested in finding the probability distribution over its possible values \( x \), denoted as \( p_X(x) \).

• It is important to distinguish between:

  • \( x \), which represents a specific value the variable can take (e.g., \( 1,2, \dots, 6 \) for a die).

  • \( X \), which is the random variable itself, generating values \(x\) according to the probability distribution \(p(x)\).

Histogramming

• Histograms provide an empirical estimate of distributions.

• Continuous distributions require density functions (PDFs), while discrete distributions use probability mass functions (PMFs).

• The bin width in histograms affects visualization, especially for discrete data.

counts, edges = np.histogram(X, range=(0,1), bins=20)

print(counts, edges)

plt.hist(X, density=True)

[3 0 0 2 3 1 2 2 1 2 5 5 4 2 2 2 1 3 7 3] [0.   0.05 0.1  0.15 0.2  0.25 0.3  0.35 0.4  0.45 0.5  0.55 0.6  0.65
 0.7  0.75 0.8  0.85 0.9  0.95 1.  ]

(array([0.60235264, 0.40156842, 1.00392106, 0.60235264, 0.60235264,
        2.00784212, 1.20470527, 0.80313685, 0.80313685, 2.00784212]),
 array([0.00192689, 0.10153631, 0.20114574, 0.30075516, 0.40036459,
        0.49997401, 0.59958344, 0.69919287, 0.79880229, 0.89841172,
        0.99802114]),
 <BarContainer object of 10 artists>)

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson

# Generate data for continuous distribution (Normal)
np.random.seed(42)
x_continuous = np.random.normal(loc=0, scale=1, size=1000)

# Generate data for discrete distribution (Poisson)
x_discrete = np.random.poisson(lam=3, size=1000)

# Define x values for theoretical curves
x_cont_range = np.linspace(-4, 4, 1000)
x_disc_range = np.arange(0, 10)

# Plot histograms
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Continuous distribution (Normal)
axes[0].hist(x_continuous, bins=30, density=True, alpha=0.6, color='b', edgecolor='black', label="Histogram")
axes[0].plot(x_cont_range, norm.pdf(x_cont_range, loc=0, scale=1), 'r-', lw=2, label="PDF")
axes[0].set_title("Continuous Distribution (Normal)")
axes[0].set_xlabel("Value")
axes[0].set_ylabel("Density")
axes[0].legend()

# Discrete distribution (Poisson)
axes[1].hist(x_discrete, bins=np.arange(11)-0.5, density=True, alpha=0.6, color='g', edgecolor='black', label="Histogram")
axes[1].scatter(x_disc_range, poisson.pmf(x_disc_range, mu=3), color='r', label="PMF", zorder=3)
axes[1].set_title("Discrete Distribution (Poisson)")
axes[1].set_xlabel("Value")
axes[1].set_ylabel("Probability")
axes[1].legend()

plt.tight_layout()
plt.show()

Expectation and Variance

• The expectation of a random variable, \( E[x] \), represents the theoretical mean, distinguishing it from the sample mean computed in simulations.

• For example, consider the difference between:

  • The average height of people computed from a sample of cities.

  • The true mean height of the entire world population.

• As the sample size increases, the sample mean converges to the expectation.

• Expectation can be applied to:

  • The variable itself (mean).

  • Any function of the variable (e.g., squared deviation for variance).

Expectation of a Random Variable

\[ E[f(x)] = \int f(x) \cdot p(x) \,dx \]

• When \( f(x) = x \), we obtain the mean, denoted by \( \mu \):

\[ E[x] = \int x \cdot p(x) \,dx = \mu \]

• Using the definition of expectation, we define variance, which quantifies the spread of \( x \).

Variance as the Expectation of Mean Fluctuations

\[ V[x] = E[(x - E[x])^2] = E[x^2] - E[x]^2 = \sigma^2 \]

• We often use the shorthand notation for variance:
  \(\sigma^2 = V[x]\), where \( \sigma \) is the standard deviation.

Binomial

• A an example of discrete distribution Binomial is defined by a Probability Mass Function (PMF)

\[P(n |p, N) = \frac{N!}{(N-n)! n!}p^n (1-p)^{N-n}\]

• \(E[n] = Np\)

• \(V[n] = 4Np(1-p)\)

Random Variable

• \(B(n, p)\) modeled by np.random.binomial(n, p, size)

Show code cell source Hide code cell source

r = np.random.binomial(n=10, p=0.6, size=2000) 

fig, ax = plt.subplots(ncols=2) 
ax[0].plot(r,  color='blue', label='trajectory')
ax[1].hist(r,  density=True, color='red',  label = 'histogram')


ax[0].set_xlabel('Samples of RN')
ax[0].set_ylabel('Values of RN')

ax[1].set_xlabel('Values of RN')
ax[1].set_ylabel('Probability Density')
fig.legend();
fig.tight_layout()

Gaussian

• A an example of continuous distribution Gaussian is defined by a Probability Distribution Function

\[P(x |\mu, \sigma) = \frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\]

• \(E[x] = \mu\)

• \(V[x] = \sigma^2\)

Random Variable

• \(N(a, b)\) modeled by np.random.normal(loc,scale, size=(N, M))

• \(N(0, 1)\) modeled by np.random.randn(N, M, P, ...)

Show code cell source Hide code cell source

# For a standard normal with sigma=1, mu=0
r = np.random.randn(200)

fig, ax = plt.subplots(ncols=2) 
ax[0].plot(r,  color='blue', label='trajectory')
ax[1].hist(r,  density=True, color='red',  label = 'histogram')


ax[0].set_xlabel('Samples of RN')
ax[0].set_ylabel('Values of RN')

ax[1].set_xlabel('Values of RN')
ax[1].set_ylabel('Probability Density')
fig.legend();
fig.tight_layout()

Uniform Distribution

• A simple example of a continuous distribution is the Uniform distribution, where all values within a given range are equally likely. It is defined by the Probability Density Function (PDF):

\[\begin{split} P(x | a, b) = \begin{cases} \frac{1}{b - a}, & a \leq x \leq b \\ 0, & \text{otherwise} \end{cases} \end{split}\]

• Expectation and Variance:

  • \(E[x] = \frac{a + b}{2}\)

  • $V[x] = \frac{(b - a)^2}{12}$

Random Variable

• \(U(a, b)\) is modeled by:

  np.random.uniform(low, high, size=(N, M))
  

• \(U(0,1)\) (standard uniform) is modeled by:

  np.random.rand(N, M, P, ...)
  

Show code cell source Hide code cell source

# For a standard uniform
r = np.random.random(200)

fig, ax = plt.subplots(ncols=2) 
ax[0].plot(r,  color='blue', label='trajectory')
ax[1].hist(r,  density=True, color='red',  label = 'histogram')


ax[0].set_xlabel('Samples of RN')
ax[0].set_ylabel('Values of RN')

ax[1].set_xlabel('Values of RN')
ax[1].set_ylabel('Probability Density')
fig.legend();
fig.tight_layout()

Exact vs sampled probability distributions

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Define range and step size
xmin, xmax, step = -4, 4, 200
dx = (xmax - xmin) / step

# Generate x values and compute normal distribution
x = np.linspace(xmin, xmax, step)
px = norm.pdf(x)  # Standard normal PDF

# Check normalization
normalization = np.sum(px * dx)
print('Normalization:', normalization)

# Generate random samples
r = np.random.randn(1000)  # Increased sample size for better histogram resolution

# Create the plot
plt.figure(figsize=(8, 5))

# Histogram of sampled data
plt.hist(r, bins=30, density=True, alpha=0.6, color='blue', edgecolor='black',
         label=f'Sampled: mean={r.mean():.2f}, var={r.var():.2f}')

# Plot theoretical normal distribution
plt.plot(x, px, 'k-', linewidth=2, label='Exact: mean=0, var=1')

# Formatting
plt.legend(loc="upper left", fontsize=10)
plt.ylabel(r'$p(x)$', fontsize=14)
plt.xlabel(r'$x$', fontsize=14)
plt.title("Comparison of Sampled Data with Normal Distribution", fontsize=12)
plt.grid(alpha=0.3)

plt.show()

Normalization: 0.9949421840184148

Learn about Transforming Random Variables

• When a random variable \(X \) is transformed by adding, multiplying by a constant, or applying a function \( Y = f(X) \), its probability distribution changes accordingly from \( p(x) \) to \( p(y) \).

• Two commonly used transformations in statistical modeling involve generating specific distributions from standard forms:

  • Generating a Gaussian (Normal) distribution from a standard normal:

    \[ N(\mu, \sigma^2) = \mu + \sigma \cdot N(0,1) \]

  • Generating a Uniform distribution from a standard uniform:

    \[ U(a, b) = (b - a) \cdot U(0,1) + a \]

• These transformations are frequently used to construct random samples from desired distributions in simulations and statistical mechanics.

Transforming Random Variables

• When transforming a random variable \( X \) to a new variable \( Y = f(X) \), the probability density functions are related by a Jacobian factor to account for how the transformation stretches or compresses the distribution:

\[ p(x) dx = p(y) dy \]

which gives:

\[ p(y) = p(x) \cdot \Bigg| \frac{dx}{dy} \Bigg| \]

• Examples of Simple Transformations:

  1. Addition: \( Y = X + a \)

     • The probability remains unchanged except for a shift:

       \[ p(y) = p(x + a) \cdot 1 \]

  2. Multiplication: \( Y = aX \)

     • The distribution scales with a factor \( \frac{1}{|a|} \):

       \[ p(y) = p(x) \cdot \frac{1}{|a|} \]

• These transformations yield useful properties:

  • Shifting the Mean:

    \[ E[X + a] = E[X] + a \]

  • Scaling the Variance:

    \[ V[aX] = a^2 V[X] \]

• Using these properties, we can generate:

  • A Gaussian (Normal) distribution from a standard normal:

    \[ N(\mu, \sigma^2) = \mu + \sigma \cdot N(0,1) \]

  • A Uniform distribution from a standard uniform:

    \[ U(a, b) = (b - a) \cdot U(0,1) + a \]

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

# Parameters
mu, sigma = 5, 2  # Mean and standard deviation for Gaussian
a, b = 2, 8       # Bounds for Uniform

# Generate standard distributions
std_normal = np.random.randn(10000)  # N(0,1)
std_uniform = np.random.rand(10000)  # U(0,1)

# Transform distributions
normal_dist = mu + sigma * std_normal  # N(mu, sigma^2)
uniform_dist = a + (b - a) * std_uniform  # U(a, b)

# Plot Distributions
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Normal Distribution
axes[0].hist(normal_dist, bins=40, density=True, alpha=0.6, color='b', edgecolor='black')
axes[0].set_title(f"Transformed Normal Distribution N({mu}, {sigma}²)")
axes[0].set_xlabel("x")
axes[0].set_ylabel("Density")

# Uniform Distribution
axes[1].hist(uniform_dist, bins=20, density=True, alpha=0.6, color='g', edgecolor='black')
axes[1].set_title(f"Transformed Uniform Distribution U({a}, {b})")
axes[1].set_xlabel("x")
axes[1].set_ylabel("Density")

plt.tight_layout()
plt.show()

Sum of Two Random Variables

• Consider the sum of two random variables, such as:

  • The sum of numbers obtained from rolling two dice.

  • The sum of two coin flips (e.g., heads = 1, tails = 0).

  • Sum of kinetic eneries of ideal gas.

\[ X = X_1 + X_2 \]

• The sum of random variables is itself a random variable!

• We want to understand how to described the properties of summed random variables as they offer a prototype of how large systems emerge froms mall components.

• Given probability distirbution of \(X_1\) and \(X_2\) how do we find probability distribution of X?

Expectation and Variance of the Sum

• Expectation is always a linear operator, which follows from the definition of expectation and the linearity of integration:

\[ E[X_1 + X_2] = E[X_1] + E[X_2] \]

• However, variance is not generally a linear operator. To see this let us write explicit formula first:

\[ V[X_1 + X_2] = E\left[(X_1 + X_2 - E[X_1 + X_2])^2\right] \]

• Defining the mean-subtracted variables: \(Y_i = X_i - E[X_i]\) we express variance of sum in terms of variances of component random variables

\[ V[X_1 + X_2] = E\left[(X_1 - E[X_1] + X_2 - E[X_2])^2\right] = E\left[(Y_1 + Y_2)^2\right] \]

• Since \( V[X_i] = E[Y_i^2] \), this simplifies to:

\[ V[X_1 + X_2] = E[Y^2_1] + V[Y^2_2] + 2E[Y_1 Y_2] = V[X_1] + V[X_2] + 2 Cov[X_1, X_2] \]

• The cross term is called Covariance which measures the degree to which two random variables vary together:

Covariance and Correlation of Two Random Variables

\[ \text{Cov}[X_1, X_2] = E[(X_1 - E[X_1])(X_2 - E[X_2])] \]

• If \(Cov > 0 \) or \(Cov < 0 \) we have positive/negative correlation and if \(Cov=0\) the variables are uncorrelated (but not necessarily independent)

To obtain a scale-independent measure, we define the correlation coefficient:

\[ \text{Corr}[X_1, X_2] = \frac{\text{Cov}[X_1, X_2]}{\sigma_{X_1} \sigma_{X_2}} \]

• In the special case where \(X_1 \) and \(X_2 \) are independent, covariane is zero and we have additivity of variances!

\[V[X_1+X_2] = V[X_1]+V[X_2]\]

• This result is fundamental in statistical mechanics, probability theory, and the sciences, as it explains why variances add for independent random variables.

Sum of \( N \) Random Variables

• Consider a sequence of independent and identically distributed (i.i.d.) random variables, \( X_1, X_2, \ldots, X_n \).

• Since they are identically distributed, each variable has a well-defined mean \( \mu \) and variance \( \sigma^2 \).

• Our goal is to understand how the sum and mean of these variables depend on the sample size \( n \).

Sample Sum and Sample Mean

\[ S_n = \sum_{i=1}^{n} X_i, \quad M_n = \frac{1}{n} \sum_{i=1}^{n} X_i \]

• \( S_n \) is the sample sum, and \( M_n \) is the sample mean.

• These quantities fluctuate with sample size \( n \), but we expect them to converge to their expectations for large \( n \)

• Because the random variables are independent, all cross terms vanish for \(i \neq j\) between mean subtracted variables \(Y_i = X_i-E[X_i]\)

Mean and Variance of the Sum of i.i.d. Random Variables

• Expectation of the Sum:

\[ E[S_n] = E\left[ \sum_{i=1}^{n} X_i \right] = \sum_{i=1}^{n} E[X_i] = n\mu \]

• Variance of the Sum:

\[ V[S_n] = E\left[ (S_n - n\mu)^2 \right] = \Bigg[\sum_{i=1}^{n} Y_i \Bigg]^2 = \sum_{i=1}^{n} \sum_{j=1}^{n}E[Y_i Y_j] = \sum_{i=1}^{n} V[X_i] = n\sigma^2 \]

Law of Large Numbers

• For the sample mean the result of summatiion of i.i.d variables implies

\[ E[M_n] = \frac{1}{n} E[S_n] = \mu \]

\[ V[M_n] = \frac{1}{n^2} V[S_n] = \frac{\sigma^2}{n} \]

• Thus, the sample mean is an unbiased estimator of \( \mu \), and its variance decreases as \( 1/n \), meaning that the estimate becomes more stable as \( n \) increases.

Law of Large Numbers (LLN)

\[ E[M_n] \to \mu \]

\[V[M_n] \to \sigma^2 / n\]

Implication:

• The sample mean provides a reliable estimate of \(\mu\) for large \(n\).

• The variance of \(M_n\) decreases as , meaning fluctuations shrink as \(1/\sqrt{n}\).

• This justifies ensemble averaging in statistical mechanics, ensuring macroscopic observables (e.g., temperature, pressure) are stable and predictable.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

# Number of trials and runs
N, runs = int(1e5), 30

# Store fractions of heads for each trial in each run
fractions = np.zeros((runs, N))

# Simulate coin tosses
for run in range(runs):
    # Generate coin tosses (0 for tails, 1 for heads)
    tosses = np.random.randint(2, size=N)
    # Calculate cumulative sum to get the number of heads up to each trial
    cum_heads = np.cumsum(tosses)
    # Calculate fraction of heads up to each trial
    fractions[run, :] = cum_heads / np.arange(1, N+1)

# Plotting
plt.figure(figsize=(14, 8))

# Plot all runs with low opacity
for run in range(runs):
    plt.plot(fractions[run, :], color='grey', alpha=0.3)

# Highlight first run
plt.semilogx(fractions[0, :], color='blue', linewidth=2, label='Highlighted Run')

# Expected value line
plt.axhline(y=0.5, color='red', linestyle='--', label='Expected Value (0.5)')
plt.xlabel('Number of Trials')
plt.ylabel('Fraction of Heads')
plt.title('Law of Large Numbers: Fraction of Heads in Coin Tossing')
plt.legend()

<matplotlib.legend.Legend at 0x7feaa79dca90>

The Central Limit Theorem (CLT)

• Central Limit Theorem asserts that the probability distribution function or PDF of sum of random variables becomes gaussian distribution with mean \(n\mu\) and \(n\sigma^2\).

• Note that CLT is based on assumption that the mean and variance, \(\mu\) and \(\sigma^2\), are finite!. Thus, CLT does not hold for certain power-law distributed random variables.

Central Limit Theorem (CLT)

• Sum of any i.i.d variables (even if they are not gaussian) leads to normally distributed random variable (the sum \(s_n\))

\[X_1 +X_2+...+X_n \rightarrow N(n\mu, n\sigma^2)\]

• The probability density function (PDF) of \(S_n\) is approaching gaussian:

\[p(s) = \frac{1}{(2\pi n\sigma^2)^{1/2}}e^{-\frac{(s-n\mu)^2}{2 n\sigma^2}}\]

• If we subtract mean and scale the sample sum by its standard deviation we will get a standard normal distribution.

\[Z_n = \frac{S_n - n\mu}{\sqrt{n}\sigma} \rightarrow N(0, 1)\]

Show code cell source Hide code cell source

from scipy.stats import norm

# Number of coin tosses in each experiment, number of experiments
N, runs    = 100, 1000  

# Simulate coin tosses: num_experiments rows, num_tosses_per_experiment columns
tosses = np.random.randint(2, size=(N, runs))

# Calculate means of each experiment
M = np.mean(tosses, axis=0)

z = ( M-M.mean() ) / np.std(M)

# Plotting the distribution of sample means
plt.figure()
plt.hist(z, density=True, bins=30)
plt.title('Distribution of Sample Means of Coin Tosses')
plt.xlabel('Sample Mean')
plt.ylabel('Density')

zs = np.linspace(z.min(), z.max(), 1000)
plt.plot(zs, norm.pdf(zs),'k', label='mean=0, var=1')
plt.legend()

<matplotlib.legend.Legend at 0x7feaa7fb5fd0>

Example of CLT applied to random walk problem

Applying the formulas to random walk model we get mean and variance for single step

\[E[X_1] = \theta \cdot 1 + (1-\theta) \cdot (-1) = 2\theta-1\]

\[V[X_1] = E[X^2_1] - E[X_1]^2 = \theta \cdot 1^2+ (1-\theta) (-1)^2 - (2\theta-1)^2 = 4 \theta(1-\theta)\]

Since steps of a random walker are independent we can compute the variance of a total displacement by multiplying mean and varaince of a single step by N

\[E[x]=N(2\theta -1)\]

\[V[x]=N\bar{\sigma^2_1} = 4N\theta (1-\theta)\]

The variance of the mean \(\bar{x} = x/N\) would then be:

\[V[\bar{x}] = \frac{4\theta (1-\theta)}{N}\]

::

Simulating a 1D unbiased random walk

• Each random walker will be modeled by a random variable \(X_i\), assuming +1 or -1 values at every step. We will run N random walkers (rows) over n steps (columns)

• We then take cumulative sum over n steps thereby summing n random variables for N walkers. This will be done via a convenient np.cumsum() method.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

def rw_1d(n, N):
    """
    Simulates a 1D symmetric random walk.

    Parameters:
    n (int): Number of steps.
    N (int): Number of walkers.

    Returns:
    np.ndarray: A (n, N) array where each column represents a walker's trajectory.
    """
    
    # Generate random steps (-1 or +1) for all walkers
    steps = np.random.choice([-1, 1], size=(n, N))
    
    # Compute cumulative sum to get displacement
    rw = np.cumsum(steps, axis=0)

    # Ensure the initial position is zero
    rw = np.vstack([np.zeros(N), rw])  # Adds a row of zeros at the start

    return rw

# Example usage: Simulate and plot a few random walks
n_steps = 1000
n_walkers = 3
rw = rw_1d(n_steps, n_walkers)

plt.plot(rw)
plt.ylabel('X (displacement)')
plt.xlabel('n (steps)')
plt.title('1D Random Walk')
plt.show()

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Simulate 1D random walk
def rw_1d(n_max, N):
    steps = np.random.choice([-1, 1], size=(n_max, N))
    return np.cumsum(steps, axis=0)

# Parameters
n_max = 1000  # Time steps
N = 1000      # Number of walkers
rw = rw_1d(n_max, N)

# Define time snapshots
time_snapshots = [10, 100, 500, 900]

# Create a multi-column subplot
fig, axes = plt.subplots(nrows=2, ncols=len(time_snapshots), figsize=(15, 6))

for i, t in enumerate(time_snapshots):
    # Plot random walk trajectories
    ax = axes[0, i]
    ax.plot(rw[:, :50], alpha=0.3)  # Show 50 trajectories for clarity
    ax.axvline(x=t, color='black', lw=2)  # Mark current time step
    ax.set_xlabel('t')
    ax.set_ylabel('X')
    ax.set_title(f'Time t={t}')

    # Histogram of positions at time t
    ax_hist = axes[1, i]
    ax_hist.hist(rw[t, :], bins=30, color='orange', density=True, alpha=0.6, label=f't={t}')
    
    # Gaussian overlay
    x = np.linspace(-100, 100, 1000)
    y = stats.norm.pdf(x, 0, np.sqrt(t))
    ax_hist.plot(x, y, color='black', lw=2, label='Normal')

    ax_hist.set_xlim([-100, 100])
    ax_hist.legend()
    ax_hist.set_title(f'$\sigma/t$ = {np.var(rw[t, :])/t:.3f}')

fig.tight_layout()
plt.show()

Mean square displacement (MSD) of a random walker

• After time n number of steps (or time t) how far has random walker moved from the origin?

\[R_n = \sum^{n-1}_{i=0}X_n\]

• We quantify this by computing Mean Square Displacement (MSD). Note that the mean is computed over N number of simulated trajectories (ensemble average). Invoking central limit theorem, or simply realizing that off diagonal terms drop off we end up with the same result as in LLN.

\[ MSD(n)= \Big\langle \big ( R_n - R_0 \big)^2 \Big \rangle \sim n \]

Show code cell source Hide code cell source

n, N = 2000, 1000
rw = rw_1d(n, N)

t = np.arange(n)

R2 = (rw[:, :] - rw[0, :])**2 # Notice we subtract initial time

msd =  np.mean(R2, axis=1)    # Notice we average over N

plt.loglog(t, np.sqrt(msd), lw=3) 

plt.loglog(t, np.sqrt(t), '--')

plt.title('Compute mean square deviation of 1D random walker',fontsize=15)
plt.xlabel('Number of steps, n',fontsize=15)
plt.ylabel(r'$MSD(n)$',fontsize=15);

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

def rw_2d(n, N):
    """
    Simulates a 2D symmetric random walk.

    Parameters:
    n (int): Number of steps.
    N (int): Number of trajectories.

    Returns:
    np.ndarray: A (n+1, N, 2) array where each trajectory is stored in the last dimension.
    """
    
    # Define possible step directions (right, up, left, down)
    steps = np.array([(1, 0), (0, 1), (-1, 0), (0, -1)])
    
    # Generate random step indices and map to step directions
    random_steps = steps[np.random.choice(4, size=(n, N))]
    
    # Prepend an initial position at (0,0) for all walkers
    rw = np.zeros((n + 1, N, 2), dtype=int)
    rw[1:] = np.cumsum(random_steps, axis=0)  # Compute displacement over time

    return rw

# Example usage: Simulate and plot first three random walkers
n_steps = 1000
n_walkers = 100
traj = rw_2d(n_steps, n_walkers)

plt.plot(traj[:, :3, 0], traj[:, :3, 1])  # Plot first three random walkers
plt.xlabel('X (displacement)')
plt.ylabel('Y (displacement)')
plt.title('2D Random Walk')
plt.show()

Show code cell source Hide code cell source

# Re-import necessary libraries due to execution state reset
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import HTML

# Parameters
num_steps = 200  # Number of steps in the random walk
step_size = 1  # Step size

# Generate random steps in 3D
np.random.seed(42)
steps = np.random.choice([-step_size, step_size], size=(num_steps, 3))  # Random steps in x, y, z
positions = np.cumsum(steps, axis=0)  # Cumulative sum to get positions

# Create figure
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim([positions[:, 0].min()-1, positions[:, 0].max()+1])
ax.set_ylim([positions[:, 1].min()-1, positions[:, 1].max()+1])
ax.set_zlim([positions[:, 2].min()-1, positions[:, 2].max()+1])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Random Walk')

# Initialize the line and point
line, = ax.plot([], [], [], 'b-', lw=2)
point, = ax.plot([], [], [], 'ro')

# Update function for animation
def update(frame):
    line.set_data(positions[:frame+1, 0], positions[:frame+1, 1])
    line.set_3d_properties(positions[:frame+1, 2])
    point.set_data([positions[frame, 0]], [positions[frame, 1]])
    point.set_3d_properties([positions[frame, 2]])
    return line, point

# Create animation
ani = animation.FuncAnimation(fig, update, frames=num_steps, interval=50, blit=False)

# Convert animation to HTML
html_anim = HTML(ani.to_jshtml())
plt.close(fig)  # Prevent duplicate display

# Display the animation
html_anim

Animation size has reached 21025445 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.

Once Loop Reflect

Show code cell source Hide code cell source

import numpy as np

import matplotlib.pyplot as plt
import matplotlib.animation as animation


def random_walk(num_steps, max_step=0.05):
    """Return a 3D random walk as (num_steps, 3) array."""
    start_pos = np.random.random(3)
    steps = np.random.uniform(-max_step, max_step, size=(num_steps, 3))
    walk = start_pos + np.cumsum(steps, axis=0)
    return walk


def update_lines(num, walks, lines):
    for line, walk in zip(lines, walks):
        line.set_data_3d(walk[:num, :].T)
    return lines


# Data: 40 random walks as (num_steps, 3) arrays
num_steps = 30
walks = [random_walk(num_steps) for index in range(40)]

# Attaching 3D axis to the figure
fig = plt.figure()
ax = fig.add_subplot(projection="3d")

# Create lines initially without data
lines = [ax.plot([], [], [])[0] for _ in walks]

# Setting the Axes properties
ax.set(xlim3d=(0, 1), xlabel='X')
ax.set(ylim3d=(0, 1), ylabel='Y')
ax.set(zlim3d=(0, 1), zlabel='Z')

# Creating the Animation object
ani = animation.FuncAnimation(
    fig, update_lines, num_steps, fargs=(walks, lines), interval=100)

plt.close()
HTML(ani.to_jshtml())

Once Loop Reflect

References

The mighty little books

• “Random Walks in Biology”, H Berg (1993)

• “Physical models of Living systems”, P Nelson (2015)

More in depth

• “Simple Brownian Diffusion: An Introduction to the Standard Theoretical Models”, D Gillespie

• “Stochastic Processes for Physicists” K Jacobs

On the applied side

• Brownian Motion: Elements of Colloid Dynamics A P Philipse (2018)

Problems

Problem 1 Binomial as generator of Gaussian and Poisson distributions

• Show that in large number limit binomial distribution tends to gaussian. Show is by expanding binomial distirbution \(logp(n)\) in power series showing that terms beyond quadratic can be ignored.

• In the limit \(N\rightarrow \infty\) but for very small values of \(p \rightarrow 0\) such that \(\lambda =pN=const\) there is another distribution that better approximates Binomial distribution: \(p(x)=\frac{\lambda^k}{k!}e^{-\lambda} \) It is known as Poisson distribution.
  Poisson distribution is an excellent approximation for probabilities of rare events. Such as, infrequently firing neurons in the brain, radioactive decay events of Plutonium or rains in the desert.
  Derive Poisson distribution by taking the limit of \(p\rightarrow 0\) in binomial distribution.

• Using numpy and matplotlib plot binomial probability distribution against Gaussian and Poisson distributions for different values of N=(10,100,1000,10000).
  

• For a value N=10000 do four plots with the following values p=0.0001, 0.001, 0.01, 0.1. You can use subplot functionality to make a pretty 4 column plot. (See plotting module)

fig, ax =  plt.subplots(nrows=1, ncols=4)
ax[0].plot()
ax[1].plot()
ax[2].plot()
ax[3].plot()

Problem-2 Confined diffusion.

Simulate 2D random walk in a circular confinement. Re-write 2D random walk code to simulate diffusion of a particle which is stuck inside a sphere. Study how root mean square deviation of position scales with time.

• Carry out simulations for different confinement sizes.

• Make plots of simulated trajectories.

Problem-3 Return to the origin!

• Simulate random walk in 1D and 2D for a different number of steps \(N=10, 10^2,10^3, 10^4, 10^5\)

• Compute average number of returns to the origin \(\langle n_{orig} \rangle\). That is number of times a random walker returns to the origin \(0\) for 1D or (0,0)$ for 2D . You may want to use some 1000 trajectories to obtain average.

• Plot how \(\langle n_{orig} \rangle\) depends on number of steps N for 1D and 2D walker.

Problem-4 Breaking the CLT; Cauchy vs Normal random walk in 2D

For this problem we are going to simulate two kinds of random walks in continuum space (not lattice): Levy flights and Normal distributd random walk.

To simulate a 2D continuum space random walk we need to generate random step sizes \(r_x\), \(r_y\). Also you will need unifrom random namber to sample angles in 2D giving you a conitnuum random walk in 2D space: \(x = r_x sin\theta\) and \(y=r_ycos\theta\)

• Normally: \(r\sim N(0,1)\)

• Cauchy distribution (long tails, infinite variance) \(r\sim Cauchy(0,1)\)

• Unform angles \(\theta \sim U(0,1)\)

Visualize random walk using matplotlib and study statistics of random walkers the way that is done for normal random walk/brownian motion examples!

Problem-5 Continuous time random walk (CTRW)

Simulate 1D random walk but instead of picking times at regular intervals pick them from exponential distribution.
Hint: you may want to use random variables from scipy.stats.exp

scipy.stats.expon

Study the root mean square deviation as a function of exponential decay parameter \(\lambda\) of exponential distribution \(e^{-\lambda x}\).